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Three dimensional (3D) numerical evolutions of static black holes with excision are presented. 
These evolutions extend to about 8000M, where M is the mass of the black hole. This degree of 
stability is achieved by using growth-rate estimates to guide the fine tuning of the parameters in a 
multi-parameter family of symmetric hyperbolic representations of the Einstein evolution equations. 
These evolutions were performed using a fixed gauge in order to separate the intrinsic stability of 
the evolution equations from the effects of stability-enhancing gauge choices. 
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Recent studies have documented the fact that 
constraint-violating instabilities are a common (if not 
universal) feature of solutions to the Einstein evolution 
equations ^, ^. Initial data with small numerical er- 
rors on some initial Cauchy surface will typically evolve 
to a solution in which the constraints grow exponen- 
tially with time. Black-hole spacetimes that are evolved 
in full 3D (without symmetry) with a fixed gauge us- 
ing one of the "standard" formulations of the evolution 
equations (e.g. ADM|, | or BSSN[| g]) have instabil- 
ities of this type that become unphysical (e.g. because 
the constraints become large) on a timescale of about 
lOOA/ ||, H, where M is the mass of the black hole. 
Several studies have shown that changing the evolution 
equations by adding multiples of the constraints and by 
changing the dynamical fields can have a significant effect 
on the growth rate of these constraint-violating instabil- 
ities [|l|, ||, ||. Such a re- formulation of the BSSN evolu- 
tion equations has allowed full 3D evolutions with fixed 
gauge to persist for about 1400Af [0. The duration of 
black hole evolutions has also been extended consider- 
ably, apparently indefinitely in some cases, by imposing 
symmetries, e.g. octant, on the solutions or by using 
an appropriate dynamical gauge ||, p^ . 

We present new results for evolving isolated static 
black holes using a multi-parameter family of symmet- 
ric hyperbolic representations of the Einstein evolution 
equations [Q. For the optimal case our evolutions ex- 
tend to about 8000M. We focus on the question of how 
the evolution equations themselves affect stability, and 
therefore we use a fixed gauge [ po| and do not impose 
any symmetries on the solutions. The fine tuning needed 
to achieve optimal stability for evolving a single black 
hole requires a special choice of the parameters in our 
representation of the evolution equations, but does not 
require any fine tuning of our numerical methods. Thus 
we expect that any numerically stable evolution code that 
solves this same system of equations with the same ini- 
tial data and boundary conditions will exhibit the same 
behavior we find here. 

We study the evolution of black-hole spacetimes using 
a particular 12-parameter family of representations of the 
Einstein evolution equations H]. This family is derived 



from the standard 3-1-1 "ADM" form of the equations by 
introducing five parameters {7, Ci Xj that densitize 
the lapse function and add multiples of the constraints 
to the evolution equations, and seven additional parame- 
ters {£, k, d, b, c, d, e} that re-define the set of dynamical 
fields. The details of the resulting evolution equations 
and the precise definitions of these various parameters 
are explained at length elsewhere ||, and will not be 
repeated here. It has been shown that a 9-parameter sub- 
family of these representations consists of strongly hyper- 
bolic evolution equations in which all of the characteristic 
speeds of the system (relative to the hypersurface- normal 
observers) have only the physical values: {0, ±1} [Q. It 
has also been shown that the evolution equations for 
an open subset of this 9-parameter family, in particular 
those representations with — | < C < 0, are symmetric 
hyperbolic Our numerical analysis here is confined 
to this 9-parameter family of symmetric hyperbolic rep- 
resentations of the Einstein evolution equations having 
physical characteristic speeds. 

Here we analyze the numerical evolution of initial 
data that represents a single isolated static black hole. 
For initial data we use a t = constant slice of the 
Schwarzschil d g eometry written in Painleve-GuUstrand 
coordinates [pJ3[ , 
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(where dVL^ is the standard metric on the unit sphere), 
plus small perturbations that are added by hand. By ex- 
plicitly inserting the same perturbations for all numerical 
resolutions, we are able to test convergence; this would 
not be the case if instead we allowed the perturbations 
to arise from machine roundoff error. The exact form of 
the perturbations is unimportant; it does not affect ei- 
ther the asymptotic growth rate of the unstable mode or 
its spatial dependence. 

We also fix the gauge for these evolutions (not just at 
the initial time but throughout the evolution) by setting 
the densitized lapse and the shift to those of Eq. (|^). 
Fixing the gauge in this way is known to be less stable 
than using a carefully selected dynamically determined 
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FIG. 1: Energy norm l|(5i5l|^/^ (per unit volume) for the most 
stable set of evolution parameters. Solid curves are 
from the full non-linear evolution code, and the dashed curve 
is from a linearized version of the code. 

gauge H, ^. However, our purpose here is to study the 
intrinsic stabihty of the evolution equations, so we choose 
to fix the gauge in this non-optimal way in order to isolate 
and emphasize this instability. 

The evolution equations are solved here using a pseu- 
dospectral collocation method (see |l^ ^ for further 
details on the implementation) on a 3D spherical shell 
extending (typically) from r = 1.9M to r = 11. 9M. This 
code utilizes the method of lines; the time integration is 
performed using a fourth-order Runge-Kutta algorithm. 
Although we use spherical coordinates, our fundamen- 
tal variables are the Cartesian components of the various 
fields. The inner boundary lies inside the event horizon; 
at this boundary all the characteristic curves are directed 
out of the domain (into the black hole) , so no boundary 
condition is required there and none is imposed ("hori- 
zon excision"). At the outer boundary we require that 
all ingoing characteristic fields be time- independent, but 
we allow all outgoing characteristic fields to propagate 
freely. 

Recent analytical work |^ has shown that the growth 
rates of the constraint-violating instabilities for the 
Painleve-GuUstrand form of the Schwarzschild geometry 
depend on just three of the nine parameters that specify 
the evolution equations, {7, C, -z}- We confine our study 
here to the dependence of this instability on the two pa- 
rameters {7, z} [ pT[ , and we fix the remaining parameters 
to the values that define System III of Ref . [|l| . 

Figure |l| shows numerical results from the evolution of 
a single black hole for the case 7 = —12, z = —0.425. 
Plotted is the energy norm (as introduced in 



FIG. 2: Solid curve shows the evolution of the integral norm 
of aU the constraints \\CkijC'''' + CkC'' + CkiijC"^"' + C^||'/^ 
(per unit volume) for the most stable set of evolution param- 
eters. Dotted curves show the individual contributions from 
the various constraints: CkijC"'^ , CkUjC'''^^ , CkC^ and C^, (in 
that order from largest to smallest at late times). 



Ref. [gj), which measures the deviation of the numerical 
solution from an exact solution that satisfies the con- 
straints. The solid curves in Figure Q represent compu- 
tations performed at different spectral resolutions (18, 
24, and 32 radial collocation points), and thus illustrate 
the convergence of our solutions. The dashed curve rep- 
resents the evolution obtained with a linearized version 
of the code, normalized so that the amplitude of the un- 
stable mode is the same as that obtained with the non- 
linear evolution. The convergence of these solutions, as 
illustrated in Fig. |], is made possible by choosing the 
same initial data, including the exact same form for the 
initial perturbation added by hand to Eq. (|^), for each 
resolution. If we had instead chosen initial data given 
by Eq. (|^) plus random perturbations (either supplied 
by numerical roundoff error or introduced by hand) we 
would not expect results using different resolutions to 
converge to the same solution. 

Figure |^ shows the evolution of the integral norm of 
the constraints (see Refs. ^ for definitions of the con- 
straint variables) for the highest-resolution case shown in 
Figure |l|. Note that at late times, most of the constraints 
in Figure || grow at the same rate (l/r « 1/275M) as 
the energy norm shown in Figure ^. The exception is 
the Hamiltonian constraint, which is much smaller than 
the other constraints, but grows at double the growth 
rate, l/r w 1/137M. Thus it appears that for the opti- 
mal choice of parameters, the unstable mode violates the 
Hamiltonian constraint only to second order in the mode 
amplitude. 
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FIG. 3; Exponential growth rates of the constraint-violating 
instabilities as a function of the parameter z (with fixed 7 = 
— 12). Points are numerically determined rates, while the solid 
curve is the approximate growth rate. 



FIG. 5: Instability growth rates as a function of the location 
of the outer boundary of the computational domain for the 
evolution parameter values 7 = — 12, z = —0.425. 
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FIG. 4: Exponential growth rates of the constraint-violating 
instabilities as a function of the parameter 7 (with fixed 2 = 
—0.425). Points are numerically determined rates, while the 
solid curve is the approximate growth rate. 



Given a numerical evolution for a particular set of pa- 
rameters, we determine the exponential growth rate by 
measuring the slope of the curve in Figs. |l] or |[ Figures ^ 
and ^ illustrate these growth rates as functions of the pa- 
rameters 7 and 2. The points in these figures represent 
numerically determined growth rates measured using the 
linearized code (which yields the same growth rates as the 
fully nonlinear code, see Figure |l| and Ref. [^). The solid 
curves represent the simple a priori estimates of these 
growth rates introduced in Ref. Q . Although the agree- 
ment between the estimates and the numerical results is 
only approximate, this agreement was good enough to 
allow us to direct our search for the most stable values 
of the parameters to the relevant region of the parame- 
ter space. The curves in these figures represent orthogo- 



nal slices of the function 1/t(7, z) through its minimum, 
l/Tmax = 1/275M, which occurs at the parameter values 
7 = —12 and z — —0.425. This minimum growth rate 
is such that constraint violations in the initial data that 
are comparable to typical machine precision (e.g. 10^^^) 
will become large (e.g. of order 0.1) when t « 10*il/. 
Figures ^ and || illustrate the full non-linear evolution 
that corresponds to this optimal choice of parameters. 

For all the cases discussed so far, the outer boundary 
radius was set at rmax = 11.9M. Figure ^ illustrates the 
dependence of the growth rate 1 /r on the location of the 
outer boundary of our computational domain, for fixed 
7 —12 and z ~ —0.425. This curve shows a sharp local 
minimum at the radius where the optimal set of evolu- 
tion parameters {7, z} was determined, strongly suggest- 
ing that these optimal values depend on the location of 
this outer boundary. We have verified this by studying in 
some detail the case where the outer boundary is located 
at Tmax = 81.9Af. There we find that the new optimal 
values of the parameters become 7 = —12 and z = —0.41, 
and the value of the growth rate at these new optimal pa- 
rameters becomes 1/r = l/333Af. This growth rate is 
about 20% smaller than that of the system whose evolu- 
tion is illustrated in Fig. ^ Thus we infer that the evo- 
lution of a single black hole in this case would extend to 
about lO^M. We also note that the optimal parameters 
for Tniax = 81. 9A/ give a value of 1/r that is about 2/3 
the value illustrated in Fig. || for this value of rmax- Con- 
siderable additional computational effort will be required 
to determine the general dependence of the optimal value 
of l/r on Tmax, and we postpone that to a future study. 
For Tijiax < 12M and for r-max > 20Af , the growth rate for 
a fixed set of evolution parameters decreases roughly like 
A/rmaxj with the constant A being about a factor of six 
larger for the case with rmax > 20Af. However, the opti- 
mal value of 1/r as a function of rmax does not scale in 
this simple way. The smallest growth rate determined in 
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our study to date is the point at 7 = —12 and z = —0.425 
with r„iax = 201. 9M, where we find 1/t = 1/570M. This 
evolution would be expected to persist for over 16,000Af . 

Finally, we note that all of the numerical evolutions 
discussed so far have placed the inner boundary of the 
domain at r,nin = 1.9M. We have also run the code with 



l.OAf and r„ 



1.5M for our best-studied case 



(7 = -12, z = -0.425, r„iax = 11.9Af) and we find that 
the growth rate is the same to three significant digits. 

In summary, we have illustrated that significant im- 
provements in the stability of numerical evolutions of 3D 
black-hole spacetimes can be achieved by a careful choice 
of the representation of the Einstein evolution equations. 
In particular we have shown that single black hole space- 
times can be evolved longer than t « SOOOAf even with 
fixed gauge. These new results also indicate that the 



outer boundary conditions may play a significant role in 
fixing the optimal formulation of the equations, as has 
been suggested by other investigations ||l^, |l^, 18, p^. 
The role of these boundary conditions will be explored 
more thoroughly in a future study. 
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